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Abstract 

We consider a model proposed by one of the authors for a type of plastic 
instability found in creep experiments which reproduces a number of exper- 
imentally observed features. The model consists of three coupled non-linear 
differential equations describing the evolution of three types of dislocations. 
The transition to the instability has been shown to be via Hopf bifurcation 
leading to limit cycle solutions with respect to physically relevant drive pa- 
rameters. Here we use reductive perturbative method to extract an amplitude 
equation of up to seventh order to obtain an approximate analytic expression 
for the order parameter. The analysis also enables us to obtain the bifurcation 
(phase) diagram of the instability. We find that while supercritical bifurca- 
tion dominates the major part of the instability region, subcritical bifurcation 
gradually takes over at one end of the region. These results are compared 
with the known experimental results. Approximate analytic expressions for the 
limit cycles for different types of bifurcations are shown to agree with their 
corresponding numerical solutions of the equations describing the model. The 
analysis also shows that high order nonlinearities are important in the prob- 
lem. This approach further allows us to map the theoretical parameters to the 
experimentally observed macroscopic quantities. 
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1 Introduction 



Instabilities in plastic flow has been an object of attention for a long time in metal- 
lurgical literature. Experimentally, there are basically three modes of deformation of 
a specimen. The best known and mostly widely studied form of the instability arises 
when the specimen is subjected to a constant rate of tensile deformation commonly 
referred to as constant strain rate test |1|, 0]. Clearly, this method of deformation 
is conceptually difficult to understand since the specimen is subjected to a prede- 
termined response ( i.e., a constant rate of deformation), and the force or the stress 
developed in the sample is sought to be measured. Under normal conditions, one 
finds a smooth stress-strain curve. However, when the system is in the regime of 
instability (i.e., for some values of the material parameters), the stress-strain curve 
exhibits repeated load drops. Each of the load drops is associated with the formation 
and propagation of dislocation bands [Q. There is another form of testing, where, the 
deformation is carried out keeping the stress rate fixed. Again, under normal condi- 
tions one finds a smooth stress-strain curve. However, when the material is deformed 
in the instability regime of the parameter space, one finds a stepped response in the 
stress-strain curve. This method of deformation is equally popular among experimen- 
talists for the study of the instability. However, conceptually the simplest form of the 
instability manifests when the material is subjected to a creep test wherein a force 
is applied and the response in the form of elongation of the specimen is measured. 
Here again, under normal conditions, the strain-time curve is smooth. Under certain 
metallurgical conditions, one sees steps on creep curve suggesting a form of instability 

1^. It is in the former two types of testing, that the plastic instability manifests 
much more easily than in the last case and hence these two modes of deformation are 
usually adopted. In contrast, the phenomenon of steps on creep curve, which is the 
subject of the present discussion, is seen in much fewer instances |[. Instabilities 
occurring in all these forms are considered to be of common origin. The phenomenon 
is referred to as the Portevin-Le Chatelier (PLC) effect or the jerky flow and is seen 
in several metals such as commercial aluminium, brass, alloys of aluminium and mag- 
nesium In the case of constant strain rate case, it is observed only in a window 
of strain rates and temperature. 

It is generally agreed that the microscopic origin of the instabilities arises due to 
the interaction of dislocations with mobile point defects and is referred to as dynamic 
strain ageing. This leads to negative strain rate characteristic of the flow stress. The 
basic idea was formulated by Cottrell ^ few decades ago. Early phenomenological 
models including Cottrell's theory and its extensions 0, || do not deal with time 
development. In contrast, techniques of dynamical systems addresses precisely this 
aspect. Recently, there has been a resurgence of interest in plastic instabilities [9- 
13] in light of introduction of new methodology borrowed from theory of dynamical 
systems. This has helped to obtain new insights hitherto not possible [10-18]. One of 
the aims of such theories is to be able to relate the microscopic dislocation mechanisms 
to the measurable macroscopic quantities. 

An attempt to understand the problem in the above perspective was made by 
Ananthakrishna and coworkers several years ago |T^, [T^. The basic idea was to de- 
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scribe the problem from the point of view of far-from-equihbrium transition, wherein 



the new temporal order could be described as a cooperative phenomenon pi] , . In 



a series of papers [ ^5] , starting from an extended Fokker-Planck equation for the 
velocity of dislocation segments, these authors arrived at a model which consisted of 
three types of dislocations and some transformations between them [|l^. The basic 
idea could be summarized by stating that limit cycle solutions arise due to nonlinear 
interaction between three different types of dislocations, suggesting a new mathemat- 
ical mechanism for the instability. Even though the spatial inhomogeneous structure 
was ignored and only the temporal oscillatory state was sought to be described, the 
model and its extensions to the case of constant strain rate test [|I3] , proved to be very 
successful in that it could explain most of the experimentally observed features such 
as the existence of bounds on strain rate for the PLC effect to occur, the negative 
strain rate sensitivity, etc. [0, |^, [1^. One other important prediction which is a di- 
rect consequence of the dynamical basis of the model is the existence of chaotic stress 
drops in a range of strain rates [p!4| , p!5|] . Recently, there has been several attempts 
which verifies this prediction [17-20]. Indeed this verification suggests that these few 
modes represent the collective degrees of freedom of dislocations. (Note that the spa- 
tially extended nature of system implies infinite degrees of freedom.) From this point 
of view, dealing with the temporal aspect appears to be justified. Description of the 
phenomenon which includes the initiation and propagation of the bands during the 
PLC effect has also been recently attempted |T6[ . 

Since the introduction of bifurcation theory into this field several years ago by our 
group |T^, [T^, several other groups have also undertaken similar lines of attack [9-12], 



25| . In the process, we feel that finer aspects of dynamical systems have been glossed 
over in this field. For instance, one often finds that casual remarks are made about 
fast and slow modes without actually going through the procedure of demonstrating 
the existence of such modes and eliminating the fast modes in favour of the slow ones 
|P5|| . In addition, under the adiabatic elimination, the resulting modes which serve 
as order parameter variables are very complicated functions of the original modes. 
Yet, hand waiving arguments have been used in building models which we believe are 
technically suspect. 

In our recent work [26| we showed how under certain conditions, one of the vari- 
ables of the model could be adiabatically eliminated since the time constant of this 
mode can be chosen to be much faster than the other two ( i.e., for low values of a 
parameter bo, see below). We then derived the equation for the order parameter for 
the reduced model. We found both supercritical and subcritical bifurcation within 
the range of applicability. We also found that the results were in good agreement 
both with the reported experimental results and with the numerical solution of the 
model. However, eliminating one of the variables entirely restricts the applicability of 
the analysis to the two-dimensional plane of the parameter space (parameters a and 
c, see below). In addition, we also found that even within the limited domain, very 
high order nonlinearities control most of the bifurcation domain. 

The purpose of the present work is to perform the analysis by keeping all the 
three modes in the model and to explore the entire instability domain spanned by 
all the three parameters (a, bo and c). In addition, this analysis should help us to 
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verify if high order nonhnearities could control part of the subcritical bifurcation as 
found in our recent analysis ^6) . This will help us to investigate the full nature of the 
bifurcation in detail. We use reductive perturbative method and extract the complex 
order parameter which is directly related to the amplitude and the frequency of the 
jumps on the creep curve. The analysis should also help us to compare the results 
with the experimental ones. The expression for the order parameter is checked by 
comparing it with the numerical solution of the model. 

In what follows (section 2) we present a brief summary of the model. In section 
3 we use the reductive perturbative method to extract the amplitude equation up to 
seventh (septic) order. This enables us not only to determine the nature of bifur- 
cation (i.e. supercritical or subcritical) exhibited by the model but also gives us an 
expression for the order parameter over most of the instability domain. In section 
4, the approximate limit cycle solution obtained through the amplitude equation is 
compared with the experimental results as well as with the numerical solution of the 
model. Section 5 contains summary and discusses our results. 



2 A Model for Steps on Creep Curve 

We start with a brief summary of the model. The details of the model can be found 



in the original references |T2[. The model consists of the mobile dislocations, the 
immobile dislocations and another type which mimics the Cottrell's type, which are 
dislocations with clouds of solute atoms. Let the corresponding densities be Nm, Nim 
and Ni, respectively. The rate equations for the densities of dislocations are: 

^ = eV^Nrr,-(3Nl-(3N^N,^ + ^Ni^-amN^, (1) 
= (3Nl-f3N,^N^-^N,,^ + aiNi, (2) 



dt 

dt 



amNm - OliNi. (3) 



The first term in Eq.(l) is the rate of production of dislocations due to cross glide 
with a rate constant 9Vm , where Vm is the velocity of the mobile dislocations which in 
general depends on some power of the applied stress, o"a. The second term refers to two 
mobile dislocations either annihilating or immobilizing. The third term also represents 
the annihilation of a mobile dislocation with an immobile one. The fourth term 
represents the remobilization of the immobile dislocations due to stress or thermal 
activation ( see jNim in Eq. 2). The last term represents the immobilization of mobile 
dislocations either due to solute atoms or due to other pinning centers, am refers to 
the concentration of the solute atoms which participate in slowing down the mobile 
dislocations. Once a mobile dislocation starts acquiring solute atoms we regard it as 
a new type of dislocation, namely the Cottrell's type iVj. This process is represented 
as an incoming term in Eq.(3). As they acquire more and more solute atoms they 
will slow down and eventually stop the dislocation entirely. At this point, they are 
considered to have transformed to Nim- This process is represented by the loss term 
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in Eq.(3) and a gain term in Eq.(2). These equations can be cast into a dimensionless 
form by using scaled variables: 

X = NU^),y = N,^{J-),z = iV,(^), r = 9V^t, (4) 
7 9Vm jam 

to get 

a; = (1 — a)x — h^x^ — xy + y, (5) 

y = bo (box^ - xy - y + azj , (6) 

z = c{x — z), (7) 

where the dot represents differentiation with respect to r while a = ^7^,60 = 
^77-, and c = Eqs.(5-7) are coupled set of nonlinear equations which support 

limit cycle solutions for a range of parameters a, Bq and c, that are physically rele- 
vant, a refers to the concentration of the solute atoms, bo refers to the reactivation of 
immobile dislocations and c to the time scales over which the slowing down occurs. 
The dependence on stress and temperature appears through Vm- We have demon- 
strated the existence of limit cycle solutions and also obtained approximate closed 



form solutions for the limit cycles II^. In addition, the model has been studied nu- 
merically. Using the Orowan equation which relates the rate of change of strain(5') to 
dislocation density and the mean velocity: S = bNmVm, with b as the Burger's vector, 
steps on the creep curve follow automatically since the densities of dislocations are 
oscillatory. Several experimental results are reproduced [O . 



3 Reductive Perturbative Approach 

We briefly outline the reductive perturbative approach to problems of formation of 
new states of order in far-from-equilibrium situations. Transitions occurring in these 
systems are quite analogous to equilibrium phase transitions. The general idea is 
to construct a 'potential like function' for the 'order parameter' like variable in the 
neighbourhood of the critical value of the drive parameter. This would permit the 
use of the methods developed in equilibrium phase transitions for further analysis. 

Near the point of Hopf bifurcation of the system (Eqs.5-7), corresponding to a 
value near the critical drive parameter, a pair of complex conjugate eigenvalues and 
another real negative eigenvalue exist for the linearized system of equations around 
the steady state. As we approach the critical value from the stable side, the real 
part of the pair of complex conjugate eigenvalues approaches zero from the negative 
side and hence the corresponding eigen directions have slow time scale. As we enter 
the instability region, these real parts become positive. In contrast, the effect of the 
change in the drive parameter on the real negative eigenvalue is negligible. Thus, 
while the two eigenvectors corresponding to the pair of complex conjugate eigenvalues 
are slow modes the eigenvector corresponding to the real negative eigenvalue is a fast 
(and decaying) mode. For this reason, the slow modes determine the formation of 
new states of order. The reductive perturbative method is a method where the slow 
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enslaving dynamics is extracted in a systematic way [27-32]. The method involves 
in first finding the critical eigenvectors corresponding to the bifurcation point and 
expressing the general solution as a linear combination of these vectors. The effect of 
the nonlinearity is handled progressively using multiple scales method. The equations 
governing the complex order parameter takes the form of Stuart-Landau equation and 
corresponds to the time- dependent Ginzburg-Landau equation for a homogeneous 
medium. (Henceforth, we will also refer to this equation as amplitude equation.) 
On the other hand, the asymptotic solution, which is a limit cycle, collapses to the 
sub-space spanned by the slow modes with no trace of the fast mode. It may be 
worth emphasizing that this method is essentially the same as reduction to center 
manifold. Indeed, the equivalence of the center manifold theory [33-35] with the 
reductive perturbation has been established |P2]. Other techniques of extracting 



amplitude equations have been devised whose end results are basically the same. For 
instance, perturbative renormalization group method |3^, ^] and its recent extension 
on the basis of envelope theory |3^ has also been developed as a tool for global 



asymptotic analysis which can be used to extract the amphtude equations. 
We start with the Eqs.(5-7). There is only one fixed point defined by: 

1 - 2a + [jl - 2af + Sbo]^' 1 
Xa = Za = TT , and ,ya = -. (8) 

Defining new variables which are deviations from the fixed point 

X = X - Xa,Y = y - ya, and , Z = z - Za, (9) 



Eqs.(5-7) become 



X = -{aX + xY + boX^ + XY), (10) 
Y = -bo(TX + 6Y -aZ -boX^ + XY), (11) 
Z = c{X-Z). (12) 



where 



a = a + 2boXa + ya-l,X = Xa-l, 

(13) 

T = ya- 2boXa,S = Xa + l. 

Eqs. (10-12) will be solved reductive perturbatively. (We note here that it is possible 
to reduce these system of equations to only two by adiabatically eliminating Eq. (10). 
This is done by noting that when 6o and c is much smaller than a, by rescaling Eq.(lO), 
it can be shown to be fast variable and hence can be adiabatically eliminated. This 
was what was done in ref |^.) Writing these equations as a matrix equation where 



the nonlinear part appears separately from the linear part, we get 



^ = Li? + iV (14) 
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where 
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Consider the stabihty of the fixed point as a function of the parameter c. The 
eigenvalues, Xi,i = 1, 0, —1, of the matrix L are determined from the cubic equation: 

A3_ta2 + PA- A = 0, (18) 

where 

T = -{a + 6bo + c), (19) 

P = 6boc + a{6bo + c) - xr6o, (20) 

and 

A = -boc[a5 + x{a-T)]. (21) 
The fixed point becomes unstable when one of the following conditions are violated 



T < 0, A < or A - PT > 0. (22) 

It is easy to show that the first two inequalities are not violated. Substituting for A, 
P and T in the third inequality of Eq. (|2^) , we get 

(a + 6bo)c^ + [{a + 5bof - xa^o]c + &o(a + 5bo){a5 - xF) < (23) 

as the condition for instability. Using the equality sign in Eq. (^) gives us the critical 
value of the drive parameter, c = cq, for a given values of a and bo. For c < cq, the 
fixed point is unstable. Since c is non-negative (negative c is unphysical), we get a 
unique Cq for the allowed pair of a and bo values within the instability. Fig.l shows a 
three-dimensional plot of the instability region involving all the three parameters of 
the model. (Note that since all the variables and the parameters are dimensionless 
quantities, we have all figures plotted in dimensionless quantities.) 

To get approximate analytical solution of Eq.(14), we follow the reductive per- 
turbative approach similar to that used by Mashiyama et al. and Richter et al. 



[pl| . We choose c = Co(l — e) with < e -C 1 and write the matrix L as a sum of two 



matrices, L = Lq + eLi, where Lq is the matrix L evaluated for c = Cq and 

(24) 
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Co 



7 



The eigenvalues of Lq are, 

Ai__i = ±«a;, and Xq = T, (25) 

where uJ^ = P [P and T being evaluated at c = cq). Taking the solution for i? as a 
growth out of the critical eigenmodes, we express it as a linear combination of these 
eigenmodes: 

-1 

R{r) = m e^^Vi + ^oe^°"ro + ^*e'*^"r^ = ^i^^'^^i (26) 

i=i 

where r^-'s are right-eigenvectors defined by Lq/^- = XjVj with r_i = r^. We also 
introduce left-eigenvectors, sj, defined by sJLq = A^sJ, where T stands for the 
transpose. Substituting this expression for R in the matrix equation, Eq. (|T^) , and 
multiplying both sides of the equation by one of the left-eigenvectors, we get an 
equation governing the corresponding amplitude: 

eA.- = e^^.^vl/.e^.- + ^ ^,,„v]>,vl>^e(^^ + A™)._ (27) 

k l,m,m<l 

Expressions for the coefficients and Qjim are given in Eq.(Cll) and Eqs.(C13,C14), 
respectively, in Appendix C. 

We express j as a power series expansion in ea: 

= e^(i) + eiljf + e^^ljf + (28) 

and introduce multiple time scales such that 

d d d n d 
dr OT OTi 0T2 

where ti = er, T2 = e^r,... Substituting these expressions for \E'j and ^ in the 
equation for the amplitudes, Eq (pT]) , we successively solve by equating terms of the 
same order in powers of e. First, terms of 0{e'^) give 

^ = 0, (30) 
implying that il)^p is constant in the time scale of r. 0[e) terms give the equation 



dr 



E ^..^^^'V^^e(^^ + ^-^^)^ (31) 

kll<k 



which upon integration gives 

e^^- = E /^.H^^'VP^e(^^ + ^')^ (32) 



k,Ll<k 
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where hjki = x^+x'-x ■ ^(^^) terms give the equation 



rj , (3) rj , (1) 

fc k,i,,i<k 

(33) 

Using compatibihty condition, we match terms that are varying on a slow time scale 
found on both sides of the equality and extract the slow dynamics: 

^ = + v^i^^mf- (34) 

Expression for r/j is given in Eq.(C12) in Appendix C. Since we are interested in 
asymptotic solutions, the equations governing the oscillatory amplitude and its 
complex conjugate \E'* are our concern. (The subscript j = l is left out from \E'i for the 
sake of brevity.) To C(e2), vEf = e2ip^^^ and, thus, Eq. (pT]) takes the form of a cubic 
Stuart-Landau equation: 

^ = e/i^ + r/|^|2^. (35) 
ar 

Note that fi{= /in) = fij. + i^ii and rji= rji) = rjr + irji are complex coefficients. ^ is 
the complex order parameter variable whose steady state solution gives its amplitude 
(squared) as 

= -e^ (36) 
and its frequency, Vt, (with = l^le**^"^) as 



VL = -^r ■ (37) 




This solution exists provided rjr is negative since /i^ is positive, r/^ is found to be 
negative over a major part of the instability region in the — a plane, as shown in 
Fig. 2 (the unshaded region). In this case, since the amplitude of the order parameter 
grows continuously in proportion to ea (see Eq.(^)), the transition is continuous ( 
'second order' type transition ) corresponding to supercritical bifurcation. There 
is a relatively small portion of the instability region, shown in the same figure in 
shades, where rj^ is found to be positive implying that the transition is discontinuous 
corresponding to subcritical bifurcation. In this regime, one has to go to quintic or 
even higher terms in the amplitude equation for obtaining an expression for the order 
parameter. We have carried out the reductive perturbative method further to derive 
the quintic 

-^^ = e/i^ + r/|^P^ + z/|^|^^ (38) 
ctr 

as well as the septic 

-^^ = e/i^ + r/|^P^ + z/|^|% + ^|^|6^ (39) 
ar 
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amplitude equations. (Expressions for the complex coefficients z/ = z/^ + zz/j and 
^ = + are given in Eqs. (A3) and (A4), respectively, in Appendix A.) The 
amplitude (squared), and frequency, fi, found from the steady state solution of 
the quintic amplitude equation are 

1*1^ 4 (-:!; ^ /(f)^ • («) 

and 

n = e^l, + r],\^^ + v,\^^ (41) 

while these same quantities found from the steady state solution of the septic ampli- 
tude equation are 

and 

n = e^^, + 7],\^^ + v,\^^ + i,\^\ (43) 

The ranges of validity of the quintic and septic amplitude equations enable us to 
describe a large portion of the subcritical bifurcation. Fig. 2 also shows the portion 
within the subcritical bifurcation where the dynamics of the system is supposed to 
be governed by the quintic (marked by dots) and septic (marked by open circles) 
amplitude equations. The rest of the domain of subcritical bifurcation where even 
higher order amplitude equation need to be considered is marked by crosses on the 
same figure. 

As can be seen from Fig. 2, there are two distinct regions of subcritical bifurcation 
located on either side of a = 0.5. The crossover from a supercritical to subcritical 
bifurcation is ffist seen for ho = 0.01 around a = 0.65 as the value of bo is reduced. 
The extent of subcritical bifurcation in the bo — a plane increases, as bo decreases. The 
validity of the amplitude equations will be examined by comparing the asymptotic 
solution derived from them with that of the limit cycle solution found by numerically 
integrating the model equations. 



4 Comparison with Experiments and Numerical 
Solutions 

4.1 Comparison with Experiments 

To start with, consider the supercritical regime where expressions for the amplitude 
and period of the limit cycle are relatively simple. Using the steady state solution of 
the cubic amplitude equation, Eqs. (|36|) and (^), the dependence of amplitude 
and period, P (cx: ^;^^), of the limit cycle on a or bo can be obtained. (Here, we will 
only consider how the two quantities, and P, depend on bo for fixed a. Similar 
analysis can be carried out for parameter a fixing bo-) The parameter 6o(= l/OVm) 
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is a function of the applied stress (da) and temperature (T). As remarked earlier, 7 
represents the stress and thermal activation, and at temperatures of interest, thermal 
activation can be completely ignored. As for the stress activation, it is a threshold 
process and therefore, can be taken to have a week dependence on stress until a 
critical value of this stress beyond which, it should show a rapid increase. One 
such functional form for 7 could be exp(— (Xc/o"), where represents the value of 
the stress beyond which the function rapidly rises. Since stress activation occurs 
at large values, one expects cXc to be large. Using the standard expression for 
Vm{(ya,T) = Vo{aa/ao)"^exp{-Em/kT) (with m > 1), we get 

&o~7M f-Ve^- (44) 

From this we see that 60 has a decreasing dependence on stress, with the term ( ^ j 
dominating up to (Tc beyond which it should increase. Clearly, bo decreases as a 
function of T. Figure 3 shows plots of and P as a function of bo. From the above 
discussion on the dependence of 60 on a and T, we see that |\E'p should increase as 
stress increases for the major part of stress value, since dc is large. On the other 
hand, it should decrease with increase in temperature. Since stress and temperature 
are measurable quantities, our predictions can be compared with experimental results. 
The amplitude and period of the limit cycle are related, respectively, to the amount 
of strain jumps and the period of the jumps on the creep curve calculated through 
the Orowan equation. There are very few experiments in this mode of testing as 
mentioned in the introduction. The only experiment where this dependence on stress 
and temperature has been measured for a limited range is that by Zagoruyko et. al. 
(in 1^). According to them, the amplitude of the strain jumps increases with stress 
while its period has a decreasing dependence on stress. Experiments from constant 
strain rate case also exhibit the same trend when the results are translated in terms 
of constant stress experiments. It is well known that the amplitude of the stress 
drops decreases with applied strain rate. In fact, even the numerical solution of the 



equations extended to the constant strain rate case predicts this behaviour . This 
implies that the amplitude of strain jumps should increase as stress increases [ 
This relation can be seen as follows. In constant strain rate case, the deformation rate 
is fixed and the stress developed in the sample is measured. When the contribution 
to the plastic strain rate increases due to increased dislocation motion (for whatever 
reasons), the stress has to fall in order to keep the applied strain rate constant. Thus, 
the relation between strain rate and stress is opposite.] Clearly, the general trend 
is consistent with the experimental results for most of the value of 60 • Zagorukuyko 
et. al. also report that the amplitude of the strain jumps increases while its period 
decreases with increase in temperature, which is consistent with our result. 

Fig. 4 shows plots of amplitude (squared) and period of the LC as a function of 
60 for a certain interval within the subcritical bifurcation. Using similar analysis as 
above, we find that while the dependence of the amplitude on bo is consistent with 
the experimental results reported by Zagoruyko et.al, the dependence of the period 
on bo is predominantly inconsistent with their report. 
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Experiments in constant strain rate case show that the stress drops are seen to 
arise both abruptly as well as continuously [0. Translating this result to the constant 
stress case, it implies that the strain jumps can arise both abruptly and continuously. 
This feature is manifest in the supercritical and the subcritical bifurcations seen in 
our calculations. 

4.2 Comparison with Numerical Solutions 

Having derived the amplitude equation, we now compare its result with the numerical 
solutions obtained via Eqs. (10-12). As stated earlier, there are two distinct types of 
solutions, namely the supercritical and the subcritical solutions. Even in the region 
of subcritical bifurcation, we have three types - first where the quintic amplitude 
equation works, second where the septic amplitude equation works and lastly, the 
rest of the instability domain where even higher order nonlinearities dominate and, 
as such, require even higher order amplitude equation. We will compare the solutions 
obtained from these with the numerically exact solutions of Eqs. (10-12). 

First consider the supercritical region. Using the steady state solution of the cubic 
amplitude equation, Eqs. (|36| , |37D , in Eq.(^), leads to an analytic expression for the 
limit cycle (LC) near bifurcation points in the domain of supercritical bifurcation. 
This is usually referred to as the secular equation. The derivation of the equations 
governing the LC are given in Appendix B. This solution can be compared with the 
numerically exact solution obtained by integrating the system of equations, Eqs. (10- 
12). Since the region of applicability in the bo — a plane is large, we choose one solution 
for large values of Bq (at the top end of the bifurcation diagram) and another value of 
bo at the lower end. Figures 5(a) and (b) show the plots of the secular equation along 
with their respective numerical solutions for two widely apart transition points within 
the domain of supercritical bifurcation. As can be observed clearly, they match very 
well. 

In a manner similar to the supercritical bifurcation, the analytic expression for the 
LC for a large portion of the domain of the subcritical bifurcation is obtained by using 
the steady state solution of the quintic or the septic amplitude equations. In this case, 
there are three distinct types of solutions: (a) first, where z/^ is negative and C,r positive 
implying that only the quintic amplitude equation has a steady state solution, (b) 
second, when z/^ and are both negative implying that both the quintic and the septic 
amplitude equations support steady state solutions, and (c) third, when z/^ is positive 
and is negative implying that the septic amplitude equation supports steady state 
solution while the quintic does not. The resulting LCs are compared with those found 
from the numerical ones. For the case (a), we have picked two different bifurcation 
points one to the right and another to the left of a = 0.5. Plots in Figs. 6(a) and 
(b) show LCs obtained from the secular equation and from the numerical ones. For 
the case (b). Fig. 7 shows three plots of LCs two of which are obtained from the 
secular equations corresponding to the quintic and septic amplitude equations while 
the third one is obtained from the numerical solution. It is clear from this figure that 
the LC derived from the septic amplitude equation is a better approximation than 
that derived from the quintic one. This shows that the 'containing' role is played not 
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only by the quintic but also by the septic (and even higher) term of the nonlinearity. 
We also find that the region over which such a situation is valid is a substantial 
portion of the validity of the quintic amplitude equation suggesting that higher order 



nonlinearities are in fact important as was indicated by our earlier work p6[. For 
the case (c), we have shown in Fig. 8 comparing the LC derived from the septic 
amplitude equation with that of the numerical one for a bifurcation point where the 
quintic amplitude equation does not hold. The agreement is reasonable except around 
the sharp turning point. 

Lastly, it is worthwhile stating that expansion of amplitude equation even up to 
septic term fails to cover the entire instability domain eventhough most of it is covered. 
This also suggests that higher order nonlinearities control the rest of the subcritical 
bifurcation domain. This feature is not usually encountered in model systems. 
The above comparison reveales the following: (i) The expression for the LC in the 
supercritical domain mimics the numerical solution very well provided e is taken to be 
small enough, (ii) In the subcritical domain also the two results generally match well 
with each other. However, for values of the parameters where the value of c is close to 
zero, (which also corresponds to small values of 60) ? higher order nonlinearities could 
supplement the contributions arising from the lower ones in determining the LC. In 
other words, if the coefficient of the next higher order has also negative real part its 
contribution may have to be included in the expression for the order parameter. 



5 Summary and Discussion 

We have carried out the reductive perturbative approach to the problem of steps 
on the creep curve regarding it as a formation of ordered state when the system is 
driven away from equilibrium. The dynamics of the system is described by two cou- 
pled amplitude equations: one for a transient order parameter \l/o and another for 
the complex order parameter ^ in the neighbourhood of the bifurcation point. The 
order parameter \1/ represents both the amplitude and the phase of the limit cycle 
solution when e > 0. Since the above derivation is valid only in the neighbourhood 
of the critical value, the expression for the order parameter \l/ is valid only for small 
e. This has been exemplified by the quite reasonable agreements between the LC 
solutions found from the analytic expression and that from the numerical integration 
of the model. We have shown that both supercritical and subcritical bifurcations are 
seen in the transitions to the instability domain. While the major part of the phase 
boundary exhibits supercritical bifurcation, subcritical bifurcation gradually dominates 
as the value of bo is reduced. The results of earlier calculation [|12| which used the 



method of relaxation oscillation and showed 'first order' type transition is consistent 
with the present one since the values of the parameters (a = .63, 60 = 10"'*) falls in 
the subcritical bifurcation domain. In our more recent work where we first adia- 



batically eliminated one of the variables and then applied the reductive perturbative 
method for the reduced model we found both supercritical and subcritical regions. 
However, due to the fact the adiabatic elimination itself was valid only for small val- 
ues of bo, the results were found to be valid in a small domain. Further, we found 
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that very high order nonhnearities were important. This result is supported by the 
present calculation. 

The present analysis shows a feature which is normally not seen in model systems, 
namely that higher order terms than quintic may control the subcritical bifurcation. 
In fact, even when the quintic term plays the 'containing role,' higher order terms 
may also contribute significantly to the 'role of containment'. Thus, in deriving 
the expression for the LC, it may be necessary to find out how higher and higher 
order nonlinear terms in the expansions in the amplitude equation contribute to the 
asymptotic solution. 

We comment here on the unusual feature of the model in the context of phase 
transitions. For conventional models, in the language of phase transitions, the 'free 
energy' is described by an expansion in power series of the order parameter up to sixth 
power. While the 'free energy' for 'second-order' phase transitions can be described 
by retaining up to fourth power in the order parameter, its is usually sufficient to 
retain the sixth power for the description of 'first-order' phase transitions (with the 
appropriate signs for the coefficients in the expansion). In the present case, however, 
we need to go to as high as eighth power (or more) of the order parameter in some 
regions of the drive parameter to describe the corresponding 'free energy' when the 
transition is of 'first-order' type. This feature is rather unusual and not found in 
conventional models of phase transitions. 

Due to the closed form expressions for the order parameter, the present calcula- 
tion helps us to map the theoretically introduced parameters to the experimentally 
measured quantities. For instance, in experiments one measures, the amplitude of the 
strain jumps on the creep curve as well as their frequency as a function of stress and 
temperature. The dependence of the amplitude of the strain jumps and its frequency 
on stress and temperature can be evaluated by using the order parameter equations in 
the Orowan's equation by noting that corresponds to the amplitude of the strain 
jumps while the frequency Q together with u gives a measure of the frequency of the 
steps. By properly relating the parameter feg to stress and temperature we found that 
for a certain range of the instability the amplitude increases while the period decreases 
as stress increases (at constant temperature) qualitatively agreeing with the reported 
experimental result of Zagorukuyko et al (see in f^). Also the amplitude increases 
while the period decreases as temperature is increased (at constant stress) which is 
again consistent with the experimental results of Zagorukuyko et. al. In addition we 
found other ranges within the instability showing various kinds of dependence of the 
order parameter on stress and temperature. 

Lastly, the present exercise has demonstrated the complicated dependence of the 
order parameter variable on the original modes. This will serve as a warning to those 
using hand waiving arguments for declaring certain modes as fast modes and others 
as slow modes in the modelling of such problems. 
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Appendix A 

In this appendix we give the expressions for the coefficients appearing in the TDGL 
equation, i.e. the expressions for n, 77, u and ^ in the equation = e/x-^ + r^l^P^-l- 
+ ^l^fj^^. Note that they are functions of the parameters a, 60 and c. They 
are given below. 

fj. = flu. (Al) 
■n= H {gijahjn^^ giakKp-r)- (A2) 

fc<a<j,7</3 

In Eq. (A2) and all the following equations summations over latin alphabets run over 
1,0, and -1 while those over greek alphabets run over ±1. Summation of Eq. (A2) 
includes only those terms which satisfy the condition a + /3 + 7 = 1. 

--i:iT^i^''+T^r\ (A3) 

k<j 

where 

T^jk^^ = m {gijaUjPfSC + giakUk/S^Sc) (A4) 
C<'5<7<|S<a 

and 

'^jk'^^ — X] (hja/stkjSC + hkaptj^sc)-) (A-5) 
/3<a,C<<5<7 

with the condition on the summations of Eqs. (A4) & (A5) to be a + /3+^ + S + (^ — 1. 
Lastly, the septic coefficient ^ is given by 

C-UT^l'^^ + Tji'^ + Tjl''^) (A6) 
k<j 



where 



and 



Tjl'^^ — ^ (gijaWj/S-ySCfXT + giakWk/3jdCp<j), (A7) 

fc<a<i,a-<p<C<5<7</3<a 

"^jk = 5Z gijkihjatSVk-iSCpa- + hkapVj^sCpcr) (A8) 

P<a,a<p<C<S<P<a 

^3,4 



^jfe ~ X] gijk{tja/3jUkSipa + tkaPiUjS(^pa) (A9) 

7<i8<a:,CT<p<C<5 

with the constraint a + /3 + 7 + 5 + C + p + o" = 1 imposed on the summations of 
Eqs. (A9), (AlO) & (All). Expressions for fijk, gjki, hjku tjkim, Ujkimn: Vjkimnp and 
Wjkimnpq are given in Appendix C. 
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Appendix B 

In this appendix, (approximate) asymptotic solution for R describing the hmit cycle 
will be derived using the steady state solution of the particular TDGL equation. 
For the instability region where cubic TDGL equation holds, expansion of ^ upto 
"0^^^ will be involved in determining R so that 

R = e4^)e^°Vo + [(e^(^) + e^(^))e^"'V+ c.c] (Bl) 

Note that the asymptotic solution does not contain 'ij;^^e^°'^ term since, Aq being 
negative, it decays in time. To 0(e2) 

^ = e5^(i) ^ ^(1) = £-5^. (B2) 

— * 

Using Eqs. (B2) and (32) in Eq. (Bl) will enable us to derive the components of R 
which are given by: 

X = I ^ I cos ((^cT + ^cc) + ^Ix I ^ TCOS (2ficr + 00.) + ^2^ I ^ I ^ + ^Sx I * TCOS (2(]cr + 9?a: ) , 

(B3) 

Y = Ay\'$\cos{ncr + 6y) + Biy\^\^cos{2Q.cT + (t)y) + B2y\'^\^ + B^y\^\^cos{2Q.cT + ipy) 

(B4) 

and 

Z = A^\^cos{n^T + e^) + Bi^\^\^cos{2n^T + 0^) + S2^|*p + B^^\^^cos{2n^T + v?J. 

(B5) 

where 

n^ = Vl + u, (B6) 
Ad = 2\ridl (B7) 



. . _i / Im{rid) \ 
Od = sm — — , 

V rid 



m 



and 



5id = 2|/imrw|, (B9) 
, . _i / Im{hiiirid)\ 

(pd = sm —- — , (BIO) 

V l^iiindl J 

B2d ^ 2\hu-ind\, (Bll) 

S3d = 2|/ii_i_irid|, (B12) 



cpd = szn — . — . (B13) 

V /ii_i_irid / 



Here the subscript d denotes x, y or z while Im denotes imaginary part of the con- 
cerned argument. Expressions for hjki and rjd are given in Appendix C. 
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For the instability region where quintic (or septic) TDGL equation holds, expansion 
of ^ upto ipj^^ (or upto i^f^) will be required in determining R so that 

N N 

R^Yl e''^'4"^e^°Vo + ^[etV'(")e*'^V+ c.c] (B14) 

n=2 n=l 

with N=4 (quintic case), N= 6 (septic case). The same procedure as done for the 
cubic case follows for the quintic as well as the septic cases. 



Appendix C 

In this appendix we first give the expressions for hkij, tkimj, Ukimnj, Vkimnpj and 
Wkimnpqj that appear as coefficients in the determination of ''pf \ '0*'^^ ^-iid 
■0^^^ respectively, such that 

e^^- = E V^i'Vr^e(^^ + ^')^ (CI) 

l<k 

E ^.■«rntfV^V£^e(^'=+^'+^-)^ (C2) 

m<l<k 

4^e'^^^ E «.fc^mntfV;'V£Vi^^e(^'=+^'+^-+^")^ (C3) 

n<Tn<l<k 

i^fe'^^ - E ^.-fe^mnptf V^'V£Vi'V^^)e(^'=+^'+^-+^"+^^)^ (C4) 

p<n<m<l<k 

and 

q<p<n<m<l<k 

Wc give only those coefficients that contribute to the asymptotic solution. Note that 
summation over latin alphabets take values 1,0 and -1 while summation over greek 
alphabets take values ±1. 

h,^, = . , . (C6) 

tja/S-y — E {djalhip-f + gjkahkfi-y) (C7) 

provided a + (3 + ^ ^ j. Otherwise, tja/s-y — 0. Cja/s-y — Xa + + \ — Aj. 

E [i9jaltl/3f5 + gjkatkPfS) + gjklihkaphi^s)] (C8) 

^ja/3jS k,l 

where Cjap^s = 

E [{QjalUlfifSC + gjkaUkp^Sc) + gjklihka/sU^SC + hiaptk-ySc)] (C9) 

^ja/i-ySC k,l 
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provided a-\-(3-\-^-\-5-\-C,^j. Otherwise, Vja/Sjsc — 0- The expression Cja/s-ysc — 
1 

jap'ySC.p X] 9jkl{hkal3UljSCp + hia/SUk-ySCp) + Qjkltka/S'y'tlSCp] (CIO) 

^jafS-ySCp k,l 

where Q = {gjaiVip-ySCp + gjkaVkp-ySCp) and Cjap^scp = x^+Xp+x^+Xs+x^+Xp-Xj ■ 
Next we give the expressions for /ijk, r)j and gjkf- 

sjrj 

i<a<fc,7</3 

/or A; 7^ /, Qjki = 'i.fxxjrkxTix + fxyjirkxTiy + rkyVi^), (C13) 

S'jfefe = fxxj^kx fxyjfkxfky-i (C'14) 



while 
where 



and 



[Si 



boSjy) 



fxyj -T-' ' (C-^^) 



r 
J 



— i^jxSjySjz), (C18) 

Sj Vj — SjxTjx + SjyVjy + Sjzfjz, (C19) 

r,x = x(co + A,), (C20) 

-(co + A,)(« + A,), (C21) 



r 



J!/ 



rjz = coX, (C22) 

= -(co + A,)(6o5 + A,), (C23) 

^jy ~ "^ix^ (C24) 

Sjz = o&oX- (C25) 
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Figure 1: The instability region determined by all the three independent parameters 
a, bo and c of the model. It is bounded by the three surfaces: co-surface (shown by 
a series of curved lines, c = plane (shown by a series of straight lines), and bo — 
plane. (Note that the Cq surface is determined from Eq. (23).) 

Figure 2: Plot of the bifurcation diagram in the a - Bq plane. The instability region 
is bounded by the parabolic-shaped curve and the bg = line. The unshaded (shaded) 
region exhibits super- (sub-) critical bifurcation. The shaded region marked by dots 
(open circles) shows the portion within the sub-critical bifurcation where quintic 
(septic) amplitude equation is supposed to hold. 

Figure 3: Plots of and P versus bo within the supercrtical domain when a — 0.45 
and e = 0.01. 

Figure 4: Plots of and P versus bo within the subcritical domain when a — 0.39 
and e = lO"! 

Figure 5: Plots of the limit cycle solutions (in X-Z plane) obtained from the secular 
equation (solid line) and that obtained from numerical integration of the model (Eqs. 
10-12) (marked with dots) when (a) a = 0.521, bo = 0.0202 and e = 0.1 and when (b) 
a = 0.6, bo = 0.01 and e = 0.001. Both bifurcation points lie within the supercritical 
bifurcation domain. 

Figure 6: Plots of the limit cycle solutions (in X-Z plane) obtained from the secular 
equation (solid line) and that obtained from numerical integration of the model (Eqs. 
10-12) (marked with dots) when (a) a = 0.647, bo = 0.009 and e = 0.01 and when 
(b) a = 0.38, bo — 0.004 and e = 0.0001. Both bifurcation points he within the 
subcritical bifurcation domain where the quintic amplitude equation holds (i.e. < 
and > 0). 

Figure 7: Plots of the limit cycle solution (in X-Z plane) obtained by using the 
secular equation derived from the quintic amplitude equation (outer solid line) , from 
the septic amplitude equation (inner solid line) and that obtained from numerical 
integration of the model (Eqs. 10-12) (marked with dots) when a = 0.375, bo = 0.004 
and e = 0.0001. The bifurcation point lies within the subcritical bifurcation domain 
where both Vr and are negative. 

Figure 8: Plots of the limit cycle solution (in X-Z plane) obtained from the secular 
equation (solid line) and that obtained from numerical integration of the model (Eqs. 
10-12) (marked with dots) when a = 0.37, bo = 0.001 and e = 0.0001. The bifurcation 
point lies within the subcritical bifurcation domain where i/^ > and < 0. 
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